from ase import Atoms, Atom
from ase.build import cut, rotate
from ase.optimize import FIRE
from ase.constraints import FixAtoms
from ase.neb import NEB
from ase.build import surface
from ase.units import J
from wizard.io import get_nth_nearest_neighbor_index, relax, dump_xyz, read_xyz, plot_band_structure
from calorine.tools import get_elastic_stiffness_tensor
from wizard.phono import PhonoCalc
from wizard.atoms import Morph
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt 
import tempfile
import numpy as np
import os
import re

unit_screw = '''135
pbc="T T T" Lattice="12.244558081272086 0.0 0.0 6.122279038869259 6.372678522968197 -0.43286219081272087 0.0 0.0 0.8658175123674912" Properties=species:S:1:pos:R:3:mass:R:1 config_type=initial screw
Fe  1.22556184e+01  9.50459364e-01  1.89626148e-03  5.58450000e+01
Fe  1.64359717e+00  9.50575972e-01  3.76147622e-01  5.58450000e+01
Fe  1.30718728e+01  9.50554770e-01  6.18908410e-01  5.58450000e+01
Fe  1.23547703e+00  2.42473852e-01  8.24627668e-02  5.58450000e+01
Fe  4.19222615e-01  2.42363604e-01  3.79492901e-01  5.58450000e+01
Fe  2.05174912e+00  2.42509894e-01  6.48035816e-01  5.58450000e+01
Fe  1.22555124e+01  2.36668905e+00  7.49138057e-02  5.58450000e+01
Fe  1.38882332e+01  2.36660071e+00  3.99517608e-01  5.58450000e+01
Fe  1.30718728e+01  2.36666078e+00  6.78845735e-01  5.58450000e+01
Fe  1.34800000e+01  1.65862898e+00  7.64057420e-02  5.58450000e+01
Fe  1.26637102e+01  1.65862191e+00  3.34309940e-01  5.58450000e+01
Fe  2.05181625e+00  1.65857597e+00  6.74296714e-01  5.58450000e+01
Fe  1.22555477e+01  3.78275618e+00  2.31009721e-01  5.58450000e+01
Fe  1.38883746e+01  3.78250883e+00  4.53508707e-01  5.58450000e+01
Fe  1.30719435e+01  3.78265018e+00  7.63425527e-01  5.58450000e+01
Fe  1.34801060e+01  3.07462191e+00  1.42893265e-01  5.58450000e+01
Fe  1.26636749e+01  3.07471731e+00  4.45990375e-01  5.58450000e+01
Fe  1.42965018e+01  3.07450177e+00  7.14111961e-01  5.58450000e+01
Fe  1.22557244e+01  5.19872792e+00  2.51131830e-01  5.58450000e+01
Fe  1.38885512e+01  5.19844523e+00  4.84710735e-01  5.58450000e+01
Fe  1.30721201e+01  5.19855124e+00  7.96581389e-01  5.58450000e+01
Fe  1.34802473e+01  4.49053004e+00  1.94706431e-01  5.58450000e+01
Fe  1.26638516e+01  4.49067138e+00  5.14793018e-01  5.58450000e+01
Fe  1.42966784e+01  4.49038869e+00  7.52560901e-01  5.58450000e+01
Fe  6.13367491e+00  2.41968551e-01 -7.54203004e-02  5.58450000e+01
Fe  7.76639576e+00  2.41769258e-01  1.71662145e-01  5.58450000e+01
Fe  6.95003534e+00  2.41855477e-01  4.74926007e-01  5.58450000e+01
Fe  1.34804594e+01  5.90646643e+00  2.12616138e-01  5.58450000e+01
Fe  1.26640283e+01  5.90657244e+00  5.24029527e-01  5.58450000e+01
Fe  1.42967845e+01  5.90639576e+00  7.71901774e-01  5.58450000e+01
Fe  6.13381625e+00  1.65797527e+00 -1.55042792e-01  5.58450000e+01
Fe  7.76639576e+00  1.65796466e+00  7.62589505e-02  5.58450000e+01
Fe  6.95010601e+00  1.65794346e+00  3.80987346e-01  5.58450000e+01
Fe  7.35826855e+00  9.49840989e-01 -1.52905311e-01  5.58450000e+01
Fe  6.54190813e+00  9.49908127e-01  1.61391014e-01  5.58450000e+01
Fe  8.17455830e+00  9.49855124e-01  4.14453901e-01  5.58450000e+01
Fe  6.13381625e+00  3.07412721e+00  4.94367223e-01  5.58450000e+01
Fe  7.76625442e+00  3.07430742e+00  8.14280442e-01  5.58450000e+01
Fe  6.95003534e+00  3.07419435e+00  1.09598594e+00  5.58450000e+01
Fe  7.35823322e+00  2.36607420e+00 -2.71657233e-01  5.58450000e+01
Fe  6.54197880e+00  2.36603534e+00  2.99054841e-02  5.58450000e+01
Fe  8.17445230e+00  2.36616254e+00  3.01732420e-01  5.58450000e+01
Fe  6.13367491e+00  4.49042403e+00  3.58512788e-01  5.58450000e+01
Fe  7.76607774e+00  4.49070671e+00  6.99352134e-01  5.58450000e+01
Fe  6.94989399e+00  4.49056537e+00  9.70948420e-01  5.58450000e+01
Fe  7.35805654e+00  3.78243816e+00  4.57861134e-01  5.58450000e+01
Fe  6.54187279e+00  3.78229682e+00  7.20426880e-01  5.58450000e+01
Fe  8.17424028e+00  3.78257951e+00  1.04352237e+00  5.58450000e+01
Fe  6.13349823e+00  5.90681979e+00  3.02666862e-01  5.58450000e+01
Fe  7.76593640e+00  5.90706714e+00  6.22674993e-01  5.58450000e+01
Fe  6.94968198e+00  5.90696113e+00  8.99827007e-01  5.58450000e+01
Fe  7.35787986e+00  5.19883392e+00  3.62214742e-01  5.58450000e+01
Fe  6.54169611e+00  5.19869258e+00  6.30099721e-01  5.58450000e+01
Fe  8.17409894e+00  5.19893993e+00  9.47973922e-01  5.58450000e+01
Fe  2.45995760e+00  9.50544170e-01  7.51326537e-02  5.58450000e+01
Fe  4.09275618e+00  9.50318021e-01  3.12602466e-01  5.58450000e+01
Fe  3.27633569e+00  9.50448763e-01  6.32844353e-01  5.58450000e+01
Fe  3.68445230e+00  2.42400353e-01  3.09619329e-02  5.58450000e+01
Fe  2.86809187e+00  2.42491873e-01  3.42968035e-01  5.58450000e+01
Fe  4.50084806e+00  2.42281979e-01  5.76312399e-01  5.58450000e+01
Fe  2.46008127e+00  2.36649117e+00  1.13614456e-01  5.58450000e+01
Fe  4.09293286e+00  2.36622968e+00  3.81563223e-01  5.58450000e+01
Fe  3.27651237e+00  2.36636396e+00  6.84788438e-01  5.58450000e+01
Fe  3.68462898e+00  1.65834276e+00  6.40547633e-02  5.58450000e+01
Fe  2.86821201e+00  1.65846996e+00  3.74119873e-01  5.58450000e+01
Fe  4.50106007e+00  1.65819788e+00  5.96184477e-01  5.58450000e+01
Fe  1.47048057e+01  3.78236749e+00  1.53433852e-01  5.58450000e+01
Fe  4.09303887e+00  3.78222615e+00  4.94050689e-01  5.58450000e+01
Fe  1.55212014e+01  3.78226148e+00  7.23773219e-01  5.58450000e+01
Fe  3.68480565e+00  3.07424735e+00  1.49101551e-01  5.58450000e+01
Fe  1.51129329e+01  3.07437102e+00  4.22396682e-01  5.58450000e+01
Fe  4.50116608e+00  3.07415548e+00  7.53401590e-01  5.58450000e+01
Fe  1.47049470e+01  5.19833922e+00  1.79783276e-01  5.58450000e+01
Fe  1.63375972e+01  5.19833922e+00  4.48577943e-01  5.58450000e+01
Fe  1.55212721e+01  5.19830389e+00  7.45491180e-01  5.58450000e+01
Fe  1.59294346e+01  4.49024735e+00  1.54882018e-01  5.58450000e+01
Fe  1.51130742e+01  4.49028269e+00  4.51678187e-01  5.58450000e+01
Fe  4.50120141e+00  4.49024735e+00  8.26672004e-01  5.58450000e+01
Fe  8.58268551e+00  2.41762898e-01 -1.19787307e-01  5.58450000e+01
Fe  1.02151943e+01  2.41936042e-01  1.98446360e-01  5.58450000e+01
Fe  9.39897527e+00  2.41814488e-01  4.66165968e-01  5.58450000e+01
Fe  1.59293993e+01  5.90639576e+00  1.70727643e-01  5.58450000e+01
Fe  1.51131449e+01  5.90636042e+00  4.69659852e-01  5.58450000e+01
Fe  1.67456537e+01  5.90650177e+00  7.40042269e-01  5.58450000e+01
Fe  8.58261484e+00  1.65804240e+00 -2.15202110e-01  5.58450000e+01
Fe  1.02150177e+01  1.65829682e+00  1.08271551e-01  5.58450000e+01
Fe  9.39883392e+00  1.65816254e+00  3.70608869e-01  5.58450000e+01
Fe  9.80703180e+00  9.50035336e-01 -1.42397562e-01  5.58450000e+01
Fe  8.99077739e+00  9.49918728e-01  1.29010187e-01  5.58450000e+01
Fe  1.06231802e+01  9.50173145e-01  4.70215951e-01  5.58450000e+01
Fe  8.58243816e+00  3.07444170e+00  5.26565767e-01  5.58450000e+01
Fe  1.02148410e+01  3.07467491e+00  7.98103982e-01  5.58450000e+01
Fe  9.39862191e+00  3.07456537e+00  1.09986179e+00  5.58450000e+01
Fe  9.80681979e+00  2.36643110e+00 -2.67562240e-01  5.58450000e+01
Fe  8.99063604e+00  2.36629682e+00  1.40761166e-02  5.58450000e+01
Fe  1.06230389e+01  2.36655124e+00  3.34378682e-01  5.58450000e+01
Fe  8.58226148e+00  4.49084806e+00  4.13707477e-01  5.58450000e+01
Fe  1.02147350e+01  4.49091873e+00  6.66272682e-01  5.58450000e+01
Fe  9.39848057e+00  4.49091873e+00  9.80840085e-01  5.58450000e+01
Fe  9.80664311e+00  3.78279152e+00  4.46918816e-01  5.58450000e+01
Fe  8.99042403e+00  3.78268551e+00  7.51887650e-01  5.58450000e+01
Fe  1.06228975e+01  3.78282686e+00  9.82516332e-01  5.58450000e+01
Fe  8.58219081e+00  5.90713781e+00  3.36304014e-01  5.58450000e+01
Fe  1.02148410e+01  5.90703180e+00  6.01833085e-01  5.58450000e+01
Fe  9.39851590e+00  5.90710247e+00  9.06466640e-01  5.58450000e+01
Fe  9.80660777e+00  5.19901060e+00  3.52835124e-01  5.58450000e+01
Fe  8.99035336e+00  5.19901060e+00  6.56312527e-01  5.58450000e+01
Fe  1.06229682e+01  5.19897527e+00  9.02978028e-01  5.58450000e+01
Fe  4.90915194e+00  9.50159011e-01 -2.44821943e-02  5.58450000e+01
Fe  5.72558304e+00  9.50021201e-01  4.95325837e-01  5.58450000e+01
Fe  5.31724382e+00  2.42112721e-01  2.49619781e-01  5.58450000e+01
Fe  4.90932862e+00  2.36611661e+00  3.17701519e-02  5.58450000e+01
Fe  5.72565371e+00  2.36604947e+00  3.74561071e-01  5.58450000e+01
Fe  5.31745583e+00  1.65807774e+00  2.19367194e-01  5.58450000e+01
Fe  4.90936396e+00  3.78219081e+00  2.71519223e-01  5.58450000e+01
Fe  5.72565371e+00  3.78222615e+00  9.46765346e-01  5.58450000e+01
Fe  5.31752650e+00  3.07410954e+00  6.29706413e-01  5.58450000e+01
Fe  1.71538516e+01  5.19840989e+00  1.53790767e-01  5.58450000e+01
Fe  5.72547703e+00  5.19851590e+00  8.86396993e-01  5.58450000e+01
Fe  5.31745583e+00  4.49028269e+00  5.93475452e-01  5.58450000e+01
Fe  1.10313781e+01  2.42071731e-01 -5.77645689e-02  5.58450000e+01
Fe  1.18475972e+01  2.42225442e-01  5.57903223e-01  5.58450000e+01
Fe  1.75618728e+01  5.90660777e+00  4.45296166e-01  5.58450000e+01
Fe  1.10312014e+01  1.65844876e+00 -1.17698110e-01  5.58450000e+01
Fe  1.18474205e+01  1.65854770e+00  5.57389357e-01  5.58450000e+01
Fe  1.14393993e+01  9.50335689e-01  2.35297954e-01  5.58450000e+01
Fe  1.10310601e+01  3.07473498e+00  4.52426558e-01  5.58450000e+01
Fe  1.18473852e+01  3.07476325e+00  7.94697597e-01  5.58450000e+01
Fe  1.14392580e+01  2.36664311e+00  2.01343431e-01  5.58450000e+01
Fe  1.10310954e+01  4.49091873e+00  3.32051074e-01  5.58450000e+01
Fe  1.18474205e+01  4.49081272e+00  8.51767343e-01  5.58450000e+01
Fe  1.14392226e+01  3.78282686e+00  6.07582721e-01  5.58450000e+01
Fe  1.10312367e+01  5.90689046e+00  2.89698996e-01  5.58450000e+01
Fe  1.18476325e+01  5.90674912e+00  8.39422936e-01  5.58450000e+01
Fe  1.14393286e+01  5.19886926e+00  5.77818795e-01  5.58450000e+01
135
pbc="T T T" Lattice="12.244558081272086 0.0 0.0 6.122279038869259 6.372678522968197 -0.43286219081272087 0.0 0.0 0.8658175123674912" Properties=species:S:1:pos:R:3:mass:R:1 config_type=initial screw
Fe  1.22556184e+01  9.50459364e-01 -5.30316396e-02  5.58450000e+01
Fe  1.64359717e+00  9.50575972e-01  3.84681272e-01  5.58450000e+01
Fe  1.30718728e+01  9.50554770e-01  5.79380431e-01  5.58450000e+01
Fe  1.23547703e+00  2.42473852e-01  9.10325230e-02  5.58450000e+01
Fe  4.19222615e-01  2.42363604e-01  3.85819724e-01  5.58450000e+01
Fe  2.05174912e+00  2.42509894e-01  6.59937491e-01  5.58450000e+01
Fe  1.22555124e+01  2.36668905e+00 -8.58908339e-02  5.58450000e+01
Fe  1.38882332e+01  2.36660071e+00  3.90453922e-01  5.58450000e+01
Fe  1.30718728e+01  2.36666078e+00  6.52428763e-01  5.58450000e+01
Fe  1.34800000e+01  1.65862898e+00  4.60003357e-02  5.58450000e+01
Fe  1.26637102e+01  1.65862191e+00  2.69222636e-01  5.58450000e+01
Fe  2.05181625e+00  1.65857597e+00  6.81542375e-01  5.58450000e+01
Fe  1.22555477e+01  3.78275618e+00  3.18508498e-01  5.58450000e+01
Fe  1.38883746e+01  3.78250883e+00  4.74816901e-01  5.58450000e+01
Fe  1.30719435e+01  3.78265018e+00  8.07995212e-01  5.58450000e+01
Fe  1.34801060e+01  3.07462191e+00  1.57389866e-01  5.58450000e+01
Fe  1.26636749e+01  3.07471731e+00  5.05516064e-01  5.58450000e+01
Fe  1.42965018e+01  3.07450177e+00  7.20243516e-01  5.58450000e+01
Fe  1.22557244e+01  5.19872792e+00  2.88990587e-01  5.58450000e+01
Fe  1.38885512e+01  5.19844523e+00  5.07890378e-01  5.58450000e+01
Fe  1.30721201e+01  5.19855124e+00  8.28192837e-01  5.58450000e+01
Fe  1.34802473e+01  4.49053004e+00  2.26083986e-01  5.58450000e+01
Fe  1.26638516e+01  4.49067138e+00  5.62919452e-01  5.58450000e+01
Fe  1.42966784e+01  4.49038869e+00  7.71912795e-01  5.58450000e+01
Fe  6.13367491e+00  2.41968551e-01 -3.78670035e-02  5.58450000e+01
Fe  7.76639576e+00  2.41769258e-01  1.87094269e-01  5.58450000e+01
Fe  6.95003534e+00  2.41855477e-01  5.02799035e-01  5.58450000e+01
Fe  1.34804594e+01  5.90646643e+00  2.35277198e-01  5.58450000e+01
Fe  1.26640283e+01  5.90657244e+00  5.50642544e-01  5.58450000e+01
Fe  1.42967845e+01  5.90639576e+00  7.89703198e-01  5.58450000e+01
Fe  6.13381625e+00  1.65797527e+00 -6.68868092e-02  5.58450000e+01
Fe  7.76639576e+00  1.65796466e+00  9.29964806e-02  5.58450000e+01
Fe  6.95010601e+00  1.65794346e+00  4.23605011e-01  5.58450000e+01
Fe  7.35826855e+00  9.49840989e-01 -1.26249014e-01  5.58450000e+01
Fe  6.54190813e+00  9.49908127e-01  2.08115102e-01  5.58450000e+01
Fe  8.17455830e+00  9.49855124e-01  4.24886724e-01  5.58450000e+01
Fe  6.13381625e+00  3.07412721e+00  3.34206424e-01  5.58450000e+01
Fe  7.76625442e+00  3.07430742e+00  8.06881671e-01  5.58450000e+01
Fe  6.95003534e+00  3.07419435e+00  1.07038617e+00  5.58450000e+01
Fe  7.35823322e+00  2.36607420e+00 -2.58275392e-01  5.58450000e+01
Fe  6.54197880e+00  2.36603534e+00  8.88180459e-02  5.58450000e+01
Fe  8.17445230e+00  2.36616254e+00  3.05591286e-01  5.58450000e+01
Fe  6.13367491e+00  4.49042403e+00  3.03373078e-01  5.58450000e+01
Fe  7.76607774e+00  4.49070671e+00  6.81375862e-01  5.58450000e+01
Fe  6.94989399e+00  4.49056537e+00  9.34402431e-01  5.58450000e+01
Fe  7.35805654e+00  3.78243816e+00  4.30641739e-01  5.58450000e+01
Fe  6.54187279e+00  3.78229682e+00  6.56047265e-01  5.58450000e+01
Fe  8.17424028e+00  3.78257951e+00  1.03434551e+00  5.58450000e+01
Fe  6.13349823e+00  5.90681979e+00  2.75731866e-01  5.58450000e+01
Fe  7.76593640e+00  5.90706714e+00  6.10461269e-01  5.58450000e+01
Fe  6.94968198e+00  5.90696113e+00  8.78990304e-01  5.58450000e+01
Fe  7.35787986e+00  5.19883392e+00  3.40503127e-01  5.58450000e+01
Fe  6.54169611e+00  5.19869258e+00  5.96665148e-01  5.58450000e+01
Fe  8.17409894e+00  5.19893993e+00  9.38643668e-01  5.58450000e+01
Fe  2.45995760e+00  9.50544170e-01  8.78905194e-02  5.58450000e+01
Fe  4.09275618e+00  9.50318021e-01  3.44986039e-01  5.58450000e+01
Fe  3.27633569e+00  9.50448763e-01  6.52858452e-01  5.58450000e+01
Fe  3.68445230e+00  2.42400353e-01  5.49554028e-02  5.58450000e+01
Fe  2.86809187e+00  2.42491873e-01  3.59841640e-01  5.58450000e+01
Fe  4.50084806e+00  2.42281979e-01  6.08985958e-01  5.58450000e+01
Fe  2.46008127e+00  2.36649117e+00  1.17018915e-01  5.58450000e+01
Fe  4.09293286e+00  2.36622968e+00  3.96564912e-01  5.58450000e+01
Fe  3.27651237e+00  2.36636396e+00  6.91135647e-01  5.58450000e+01
Fe  3.68462898e+00  1.65834276e+00  8.60592473e-02  5.58450000e+01
Fe  2.86821201e+00  1.65846996e+00  3.86077064e-01  5.58450000e+01
Fe  4.50106007e+00  1.65819788e+00  6.42120208e-01  5.58450000e+01
Fe  1.47048057e+01  3.78236749e+00  1.64990251e-01  5.58450000e+01
Fe  4.09303887e+00  3.78222615e+00  4.62694431e-01  5.58450000e+01
Fe  1.55212014e+01  3.78226148e+00  7.30760064e-01  5.58450000e+01
Fe  3.68480565e+00  3.07424735e+00  1.39764678e-01  5.58450000e+01
Fe  1.51129329e+01  3.07437102e+00  4.25679519e-01  5.58450000e+01
Fe  4.50116608e+00  3.07415548e+00  7.26136731e-01  5.58450000e+01
Fe  1.47049470e+01  5.19833922e+00  1.96053622e-01  5.58450000e+01
Fe  1.63375972e+01  5.19833922e+00  4.56803590e-01  5.58450000e+01
Fe  1.55212721e+01  5.19830389e+00  7.56942569e-01  5.58450000e+01
Fe  1.59294346e+01  4.49024735e+00  1.63095042e-01  5.58450000e+01
Fe  1.51130742e+01  4.49028269e+00  4.63986519e-01  5.58450000e+01
Fe  4.50120141e+00  4.49024735e+00  7.85898364e-01  5.58450000e+01
Fe  8.58268551e+00  2.41762898e-01 -1.16396466e-01  5.58450000e+01
Fe  1.02151943e+01  2.41936042e-01  1.77832883e-01  5.58450000e+01
Fe  9.39897527e+00  2.41814488e-01  4.57804336e-01  5.58450000e+01
Fe  1.59293993e+01  5.90639576e+00  1.80902000e-01  5.58450000e+01
Fe  1.51131449e+01  5.90636042e+00  4.83163912e-01  5.58450000e+01
Fe  1.67456537e+01  5.90650177e+00  7.47763647e-01  5.58450000e+01
Fe  8.58261484e+00  1.65804240e+00 -2.12096375e-01  5.58450000e+01
Fe  1.02150177e+01  1.65829682e+00  8.20537527e-02  5.58450000e+01
Fe  9.39883392e+00  1.65816254e+00  3.62093770e-01  5.58450000e+01
Fe  9.80703180e+00  9.50035336e-01 -1.59393035e-01  5.58450000e+01
Fe  8.99077739e+00  9.49918728e-01  1.26181739e-01  5.58450000e+01
Fe  1.06231802e+01  9.50173145e-01  4.34981247e-01  5.58450000e+01
Fe  8.58243816e+00  3.07444170e+00  5.25298915e-01  5.58450000e+01
Fe  1.02148410e+01  3.07467491e+00  8.10979375e-01  5.58450000e+01
Fe  9.39862191e+00  3.07456537e+00  1.10345511e+00  5.58450000e+01
Fe  9.80681979e+00  2.36643110e+00 -2.74654629e-01  5.58450000e+01
Fe  8.99063604e+00  2.36629682e+00  1.30344170e-02  5.58450000e+01
Fe  1.06230389e+01  2.36655124e+00  3.09569986e-01  5.58450000e+01
Fe  8.58226148e+00  4.49084806e+00  4.10036442e-01  5.58450000e+01
Fe  1.02147350e+01  4.49091873e+00  6.91812908e-01  5.58450000e+01
Fe  9.39848057e+00  4.49091873e+00  9.90385618e-01  5.58450000e+01
Fe  9.80664311e+00  3.78279152e+00  4.62861095e-01  5.58450000e+01
Fe  8.99042403e+00  3.78268551e+00  7.54380364e-01  5.58450000e+01
Fe  1.06228975e+01  3.78282686e+00  1.02376842e+00  5.58450000e+01
Fe  8.58219081e+00  5.90713781e+00  3.33462382e-01  5.58450000e+01
Fe  1.02148410e+01  5.90703180e+00  6.17506403e-01  5.58450000e+01
Fe  9.39851590e+00  5.90710247e+00  9.13051212e-01  5.58450000e+01
Fe  9.80660777e+00  5.19901060e+00  3.67251611e-01  5.58450000e+01
Fe  8.99035336e+00  5.19901060e+00  6.58764286e-01  5.58450000e+01
Fe  1.06229682e+01  5.19897527e+00  9.29686053e-01  5.58450000e+01
Fe  4.90915194e+00  9.50159011e-01  2.50834558e-02  5.58450000e+01
Fe  5.72558304e+00  9.50021201e-01  5.54160435e-01  5.58450000e+01
Fe  5.31724382e+00  2.42112721e-01  2.88735739e-01  5.58450000e+01
Fe  4.90932862e+00  2.36611661e+00  9.36230495e-02  5.58450000e+01
Fe  5.72565371e+00  2.36604947e+00  6.12050307e-01  5.58450000e+01
Fe  5.31745583e+00  1.65807774e+00  3.09149527e-01  5.58450000e+01
Fe  4.90936396e+00  3.78219081e+00  2.04510809e-01  5.58450000e+01
Fe  5.72565371e+00  3.78222615e+00  8.46701784e-01  5.58450000e+01
Fe  5.31752650e+00  3.07410954e+00  4.63755746e-01  5.58450000e+01
Fe  1.71538516e+01  5.19840989e+00  1.59847883e-01  5.58450000e+01
Fe  5.72547703e+00  5.19851590e+00  8.46819986e-01  5.58450000e+01
Fe  5.31745583e+00  4.49028269e+00  5.36888230e-01  5.58450000e+01
Fe  1.10313781e+01  2.42071731e-01 -8.99341201e-02  5.58450000e+01
Fe  1.18475972e+01  2.42225442e-01  5.19658223e-01  5.58450000e+01
Fe  1.75618728e+01  5.90660777e+00  4.51234021e-01  5.58450000e+01
Fe  1.10312014e+01  1.65844876e+00 -1.80205972e-01  5.58450000e+01
Fe  1.18474205e+01  1.65854770e+00  4.59880371e-01  5.58450000e+01
Fe  1.14393993e+01  9.50335689e-01  1.81826710e-01  5.58450000e+01
Fe  1.10310601e+01  3.07473498e+00  5.09242502e-01  5.58450000e+01
Fe  1.18473852e+01  3.07476325e+00  1.02899617e+00  5.58450000e+01
Fe  1.14392580e+01  2.36664311e+00  4.58993604e-02  5.58450000e+01
Fe  1.10310954e+01  4.49091873e+00  3.77303523e-01  5.58450000e+01
Fe  1.18474205e+01  4.49081272e+00  9.08940746e-01  5.58450000e+01
Fe  1.14392226e+01  3.78282686e+00  6.93504258e-01  5.58450000e+01
Fe  1.10312367e+01  5.90689046e+00  3.12953975e-01  5.58450000e+01
Fe  1.18476325e+01  5.90674912e+00  8.66693240e-01  5.58450000e+01
Fe  1.14393286e+01  5.19886926e+00  6.14084770e-01  5.58450000e+01
135
pbc="T T T" Lattice="12.244558081272086 0.0 0.0 6.122279038869259 6.372678522968197 -0.43286219081272087 0.0 0.0 0.8658175123674912" Properties=species:S:1:pos:R:3:mass:R:1 config_type=initial screw
Fe  1.22556184e+01  9.50459364e-01  6.34496466e-03  5.58450000e+01
Fe  1.64359717e+00  9.50575972e-01  3.86529611e-01  5.58450000e+01
Fe  1.30718728e+01  9.50554770e-01  6.22473092e-01  5.58450000e+01
Fe  1.23547703e+00  2.42473852e-01  9.33693039e-02  5.58450000e+01
Fe  4.19222615e-01  2.42363604e-01  3.87849014e-01  5.58450000e+01
Fe  2.05174912e+00  2.42509894e-01  6.62655753e-01  5.58450000e+01
Fe  1.22555124e+01  2.36668905e+00  7.58668057e-02  5.58450000e+01
Fe  1.38882332e+01  2.36660071e+00  4.00127710e-01  5.58450000e+01
Fe  1.30718728e+01  2.36666078e+00  6.79598795e-01  5.58450000e+01
Fe  1.34800000e+01  1.65862898e+00  7.83994452e-02  5.58450000e+01
Fe  1.26637102e+01  1.65862191e+00  3.36790042e-01  5.58450000e+01
Fe  2.05181625e+00  1.65857597e+00  6.82767806e-01  5.58450000e+01
Fe  1.22555477e+01  3.78275618e+00  2.28218152e-01  5.58450000e+01
Fe  1.38883746e+01  3.78250883e+00  4.51707463e-01  5.58450000e+01
Fe  1.30719435e+01  3.78265018e+00  7.61209269e-01  5.58450000e+01
Fe  1.34801060e+01  3.07462191e+00  1.42218173e-01  5.58450000e+01
Fe  1.26636749e+01  3.07471731e+00  4.45146873e-01  5.58450000e+01
Fe  1.42965018e+01  3.07450177e+00  7.13559505e-01  5.58450000e+01
Fe  1.22557244e+01  5.19872792e+00  2.45291110e-01  5.58450000e+01
Fe  1.38885512e+01  5.19844523e+00  4.80797792e-01  5.58450000e+01
Fe  1.30721201e+01  5.19855124e+00  7.91843004e-01  5.58450000e+01
Fe  1.34802473e+01  4.49053004e+00  1.91490700e-01  5.58450000e+01
Fe  1.26638516e+01  4.49067138e+00  5.10823329e-01  5.58450000e+01
Fe  1.42966784e+01  4.49038869e+00  7.49905601e-01  5.58450000e+01
Fe  6.13367491e+00  2.41968551e-01 -3.07034382e-02  5.58450000e+01
Fe  7.76639576e+00  2.41769258e-01  1.99292576e-01  5.58450000e+01
Fe  6.95003534e+00  2.41855477e-01  5.12026113e-01  5.58450000e+01
Fe  1.34804594e+01  5.90646643e+00  2.07407512e-01  5.58450000e+01
Fe  1.26640283e+01  5.90657244e+00  5.17747898e-01  5.58450000e+01
Fe  1.42967845e+01  5.90639576e+00  7.67523986e-01  5.58450000e+01
Fe  6.13381625e+00  1.65797527e+00 -6.33555936e-02  5.58450000e+01
Fe  7.76639576e+00  1.65796466e+00  9.97222686e-02  5.58450000e+01
Fe  6.95010601e+00  1.65794346e+00  4.28359558e-01  5.58450000e+01
Fe  7.35826855e+00  9.49840989e-01 -1.17709226e-01  5.58450000e+01
Fe  6.54190813e+00  9.49908127e-01  2.14466417e-01  5.58450000e+01
Fe  8.17455830e+00  9.49855124e-01  4.36862643e-01  5.58450000e+01
Fe  6.13381625e+00  3.07412721e+00  3.32996403e-01  5.58450000e+01
Fe  7.76625442e+00  3.07430742e+00  8.04508788e-01  5.58450000e+01
Fe  6.95003534e+00  3.07419435e+00  1.06873894e+00  5.58450000e+01
Fe  7.35823322e+00  2.36607420e+00 -2.56307131e-01  5.58450000e+01
Fe  6.54197880e+00  2.36603534e+00  9.02274205e-02  5.58450000e+01
Fe  8.17445230e+00  2.36616254e+00  3.08531148e-01  5.58450000e+01
Fe  6.13367491e+00  4.49042403e+00  2.97825643e-01  5.58450000e+01
Fe  7.76607774e+00  4.49070671e+00  6.71327004e-01  5.58450000e+01
Fe  6.94989399e+00  4.49056537e+00  9.27077208e-01  5.58450000e+01
Fe  7.35805654e+00  3.78243816e+00  4.25033580e-01  5.58450000e+01
Fe  6.54187279e+00  3.78229682e+00  6.51976770e-01  5.58450000e+01
Fe  8.17424028e+00  3.78257951e+00  1.02616188e+00  5.58450000e+01
Fe  6.13349823e+00  5.90681979e+00  2.67391258e-01  5.58450000e+01
Fe  7.76593640e+00  5.90706714e+00  5.97144004e-01  5.58450000e+01
Fe  6.94968198e+00  5.90696113e+00  8.68532127e-01  5.58450000e+01
Fe  7.35787986e+00  5.19883392e+00  3.29933064e-01  5.58450000e+01
Fe  6.54169611e+00  5.19869258e+00  5.88563039e-01  5.58450000e+01
Fe  8.17409894e+00  5.19893993e+00  9.24481894e-01  5.58450000e+01
Fe  2.45995760e+00  9.50544170e-01  9.00631095e-02  5.58450000e+01
Fe  4.09275618e+00  9.50318021e-01  3.48121399e-01  5.58450000e+01
Fe  3.27633569e+00  9.50448763e-01  6.55447452e-01  5.58450000e+01
Fe  3.68445230e+00  2.42400353e-01  5.87700283e-02  5.58450000e+01
Fe  2.86809187e+00  2.42491873e-01  3.63040449e-01  5.58450000e+01
Fe  4.50084806e+00  2.42281979e-01  6.13605173e-01  5.58450000e+01
Fe  2.46008127e+00  2.36649117e+00  1.17468675e-01  5.58450000e+01
Fe  4.09293286e+00  2.36622968e+00  3.97224413e-01  5.58450000e+01
Fe  3.27651237e+00  2.36636396e+00  6.91675325e-01  5.58450000e+01
Fe  3.68462898e+00  1.65834276e+00  8.78160671e-02  5.58450000e+01
Fe  2.86821201e+00  1.65846996e+00  3.87532633e-01  5.58450000e+01
Fe  4.50106007e+00  1.65819788e+00  6.44281767e-01  5.58450000e+01
Fe  1.47048057e+01  3.78236749e+00  1.51941519e-01  5.58450000e+01
Fe  4.09303887e+00  3.78222615e+00  4.60752873e-01  5.58450000e+01
Fe  1.55212014e+01  3.78226148e+00  7.22516852e-01  5.58450000e+01
Fe  3.68480565e+00  3.07424735e+00  1.39171495e-01  5.58450000e+01
Fe  1.51129329e+01  3.07437102e+00  4.21936286e-01  5.58450000e+01
Fe  4.50116608e+00  3.07415548e+00  7.25404629e-01  5.58450000e+01
Fe  1.47049470e+01  5.19833922e+00  1.76502113e-01  5.58450000e+01
Fe  1.63375972e+01  5.19833922e+00  4.46181078e-01  5.58450000e+01
Fe  1.55212721e+01  5.19830389e+00  7.42702975e-01  5.58450000e+01
Fe  1.59294346e+01  4.49024735e+00  1.52986406e-01  5.58450000e+01
Fe  1.51130742e+01  4.49028269e+00  4.49450014e-01  5.58450000e+01
Fe  4.50120141e+00  4.49024735e+00  7.82425792e-01  5.58450000e+01
Fe  8.58268551e+00  2.41762898e-01 -9.98429435e-02  5.58450000e+01
Fe  1.02151943e+01  2.41936042e-01  2.09294424e-01  5.58450000e+01
Fe  9.39897527e+00  2.41814488e-01  4.80699502e-01  5.58450000e+01
Fe  1.59293993e+01  5.90639576e+00  1.67524806e-01  5.58450000e+01
Fe  1.51131449e+01  5.90636042e+00  4.65935565e-01  5.58450000e+01
Fe  1.67456537e+01  5.90650177e+00  7.37261191e-01  5.58450000e+01
Fe  8.58261484e+00  1.65804240e+00 -2.01919972e-01  5.58450000e+01
Fe  1.02150177e+01  1.65829682e+00  1.14031307e-01  5.58450000e+01
Fe  9.39883392e+00  1.65816254e+00  3.79014155e-01  5.58450000e+01
Fe  9.80703180e+00  9.50035336e-01 -1.32082060e-01  5.58450000e+01
Fe  8.99077739e+00  9.49918728e-01  1.43829505e-01  5.58450000e+01
Fe  1.06231802e+01  9.50173145e-01  4.77735703e-01  5.58450000e+01
Fe  8.58243816e+00  3.07444170e+00  5.21589025e-01  5.58450000e+01
Fe  1.02148410e+01  3.07467491e+00  7.96089403e-01  5.58450000e+01
Fe  9.39862191e+00  3.07456537e+00  1.09685244e+00  5.58450000e+01
Fe  9.80681979e+00  2.36643110e+00 -2.65120456e-01  5.58450000e+01
Fe  8.99063604e+00  2.36629682e+00  1.78944346e-02  5.58450000e+01
Fe  1.06230389e+01  2.36655124e+00  3.36073558e-01  5.58450000e+01
Fe  8.58226148e+00  4.49084806e+00  3.95600198e-01  5.58450000e+01
Fe  1.02147350e+01  4.49091873e+00  6.57512707e-01  5.58450000e+01
Fe  9.39848057e+00  4.49091873e+00  9.68554145e-01  5.58450000e+01
Fe  9.80664311e+00  3.78279152e+00  4.40021693e-01  5.58450000e+01
Fe  8.99042403e+00  3.78268551e+00  7.41450148e-01  5.58450000e+01
Fe  1.06228975e+01  3.78282686e+00  9.77641099e-01  5.58450000e+01
Fe  8.58219081e+00  5.90713781e+00  3.16314194e-01  5.58450000e+01
Fe  1.02148410e+01  5.90703180e+00  5.89744053e-01  5.58450000e+01
Fe  9.39851590e+00  5.90710247e+00  8.90972936e-01  5.58450000e+01
Fe  9.80660777e+00  5.19901060e+00  3.40320184e-01  5.58450000e+01
Fe  8.99035336e+00  5.19901060e+00  6.39328293e-01  5.58450000e+01
Fe  1.06229682e+01  5.19897527e+00  8.93511940e-01  5.58450000e+01
Fe  4.90915194e+00  9.50159011e-01  2.89538693e-02  5.58450000e+01
Fe  5.72558304e+00  9.50021201e-01  5.59049548e-01  5.58450000e+01
Fe  5.31724382e+00  2.42112721e-01  2.94429392e-01  5.58450000e+01
Fe  4.90932862e+00  2.36611661e+00  9.44470848e-02  5.58450000e+01
Fe  5.72565371e+00  2.36604947e+00  6.13108883e-01  5.58450000e+01
Fe  5.31745583e+00  1.65807774e+00  3.11872155e-01  5.58450000e+01
Fe  4.90936396e+00  3.78219081e+00  2.02095661e-01  5.58450000e+01
Fe  5.72565371e+00  3.78222615e+00  8.43617746e-01  5.58450000e+01
Fe  5.31752650e+00  3.07410954e+00  4.62829332e-01  5.58450000e+01
Fe  1.71538516e+01  5.19840989e+00  1.51709339e-01  5.58450000e+01
Fe  5.72547703e+00  5.19851590e+00  8.40453145e-01  5.58450000e+01
Fe  5.31745583e+00  4.49028269e+00  5.32554608e-01  5.58450000e+01
Fe  1.10313781e+01  2.42071731e-01 -4.94491343e-02  5.58450000e+01
Fe  1.18475972e+01  2.42225442e-01  5.64437852e-01  5.58450000e+01
Fe  1.75618728e+01  5.90660777e+00  4.42860463e-01  5.58450000e+01
Fe  1.10312014e+01  1.65844876e+00 -1.13517760e-01  5.58450000e+01
Fe  1.18474205e+01  1.65854770e+00  5.60556484e-01  5.58450000e+01
Fe  1.14393993e+01  9.50335689e-01  2.40992654e-01  5.58450000e+01
Fe  1.10310601e+01  3.07473498e+00  4.50984042e-01  5.58450000e+01
Fe  1.18473852e+01  3.07476325e+00  7.93614053e-01  5.58450000e+01
Fe  1.14392580e+01  2.36664311e+00  2.02588300e-01  5.58450000e+01
Fe  1.10310954e+01  4.49091873e+00  3.25536300e-01  5.58450000e+01
Fe  1.18474205e+01  4.49081272e+00  8.46752640e-01  5.58450000e+01
Fe  1.14392226e+01  3.78282686e+00  6.03962074e-01  5.58450000e+01
Fe  1.10312367e+01  5.90689046e+00  2.80134827e-01  5.58450000e+01
Fe  1.18476325e+01  5.90674912e+00  8.31733530e-01  5.58450000e+01
Fe  1.14393286e+01  5.19886926e+00  5.70470011e-01  5.58450000e+01
'''

class MaterialCalculator():
    def __init__(self, atoms, calculator, symbol_info):
        atoms.calc = calculator
        relax(atoms)
        atom_energy = atoms.get_potential_energy() / len(atoms)
        self.atom_energy = atom_energy
        self.atoms = atoms
        self.calc = calculator  
        self.formula = symbol_info.formula
        self.lc = symbol_info.lattice_constant
        self.structure = symbol_info.structure
    
    def isolate_atom_energy(self):
        symbol = self.atoms.get_chemical_symbols()[0]
        atoms = Atoms(symbols = [symbol], positions=[[0,0,0]])
        atoms.calc = self.calc
        iso_atom_energy = atoms.get_potential_energy()
        with open('MaterialProperties.out', "a") as f:
            f.write(f" {self.formula:<7}Atom_Energy: {iso_atom_energy:.4f} eV\n")
        return iso_atom_energy

    def get_potential_energy(self):
        energy = self.atom_energy
        return energy

    def lattice_constant(self):
        atoms = self.atoms
        atom_energy = self.atom_energy
        cell_lengths = atoms.cell.cellpar()
        dump_xyz('MaterialProperties.xyz', atoms)
        
        output = ""
        if self.structure == 'hcp':
            output += f" {self.formula:<10}Lattice_Constants: a: {cell_lengths[0]:.4f} A    c: {cell_lengths[2]:.4f} A\n"
            output += f"{'':<11}Ground_State_Energy: {atom_energy:.4f} eV\n"
        else:
            output += f" {self.formula:<10}Lattice_Constants: {round(sum(cell_lengths[:3])/3, 3):.4f} A\n"
            output += f"{'':<11}Ground_State_Energy: {atom_energy:.4f} eV\n"
        
        with open('MaterialProperties.out', "a") as f:
            f.write(output)
        
        return output
    
    def elastic_constant(self, epsilon = 0.01):
        atoms = self.atoms.copy()
        atoms.calc = self.calc
        Cij = get_elastic_stiffness_tensor(atoms, epsilon=epsilon)
        dump_xyz('MaterialProperties.xyz', atoms)
        
        output = ""
        output += f" {self.formula:<10}C11: {Cij[0][0]:>7.2f} GPa\n"
        output += f"{'':<11}C12: {Cij[0][1]:>7.2f} GPa\n"
        output += f"{'':<11}C13: {Cij[0][2]:>7.2f} GPa\n"
        output += f"{'':<11}C33: {Cij[2][2]:>7.2f} GPa\n"
        output += f"{'':<11}C44: {Cij[3][3]:>7.2f} GPa\n"
        output += f"{'':<11}C66: {Cij[5][5]:>7.2f} GPa\n"
        
        with open('MaterialProperties.out', 'a') as f:
            f.write(output)
        
        return output
    
    def eos_curve(self):
        atoms = self.atoms.copy()
        atoms.calc = self.calc
        os.makedirs('eos_curve_out', exist_ok=True)
        os.makedirs('eos_curve_png', exist_ok=True)
        volumes, energies = [], []
        origin_cell = atoms.cell.copy()
        for scale in np.arange(0.9, 1.10, 0.01):
            atoms.set_cell(scale * origin_cell, scale_atoms = True)
            volumes.append(atoms.get_volume() / len(atoms))
            energies.append(atoms.get_potential_energy() / len(atoms))
            dump_xyz('MaterialProperties.xyz', atoms)
            
        fig, ax = plt.subplots()
        plt.subplots_adjust(left=0.15, right=0.95, top=0.95, bottom=0.15)
        font_size = 12
        ax.plot(volumes, energies, '-o')
        ax.set_xlabel('Volume(A$^3$/atom)', fontsize=font_size)
        ax.set_ylabel('Energy (eV/atom)', fontsize=font_size)
        ax.set_title(f'{self.formula} {self.structure} EOS Curve', fontsize=font_size)
        fig_path = os.path.join('eos_curve_png',f'{self.formula}_eos_curve.png')
        fig.savefig(fig_path)
        plt.close(fig)
        with open(os.path.join('eos_curve_out',f'{self.formula}_eos_curve.out'), 'w') as f:
            f.write("Volume(A^3/atom)   Energy(eV/atom)\n")
            for volume, energy in zip(volumes, energies):
                f.write(f"{volume:.2f}   {energy:.4f}\n")

        return fig_path

    def phonon_dispersion(self, special_points = None, labels_path = None):
        atoms = self.atoms.copy()
        calc = self.calc
        PhonoCalc(atoms, calc).get_band_structure(special_points=special_points, labels_path=labels_path)
        fig_path = plot_band_structure(atoms, self.formula, self.structure)
        return fig_path
             
    def formation_energy_surface(self, hkl = (1, 0, 0), layers = 10, relax_params = {}):
        atoms = self.atoms.copy()
        atom_energy = self.atom_energy
        slab = surface(atoms, hkl, layers = layers, vacuum=10) 
        Morph(slab).shuffle_symbols()
        slab.calc = self.calc
        relax(slab, **relax_params)
        box = slab.get_cell()
        S = (box[0][0] * box[1][1] - box[0][1] * box[1][0])
        slab_energy = slab.get_potential_energy()
        formation_energy = (slab_energy - atom_energy * len(slab)) / S / 2

        if self.structure == 'hcp':
            hk, l = hkl[:2], hkl[2]
            hkl_str = '-'.join(map(str, sorted(hk, reverse=True)))
            hkl_str += f'-{l}'
        else:
            hkl_str = '-'.join(map(str, sorted(hkl, reverse=True)))

        dump_xyz('MaterialProperties.xyz', slab)
        with open('MaterialProperties.out', 'a') as f:
            print(f' {self.formula:<7}{hkl_str} Surface_Energy: {formation_energy / J / 1e-20 :.4f} J/m^2', file=f)
        return formation_energy * 1000

    def formation_energy_vacancy(self, index = 0, relax_params = {}):
        atoms = self.atoms.copy()
        atoms.calc = self.calc
        atom_energy = self.atom_energy
        Morph(atoms).create_vacancy(index = index)
        relax(atoms, **relax_params)
        formation_energy = atoms.get_potential_energy() - atom_energy * len(atoms)

        dump_xyz('MaterialProperties.xyz', atoms)
        with open('MaterialProperties.out', 'a') as f:
            print(f' {self.formula:<7}Formation_Energy_Vacancy: {formation_energy:.4f} eV', file=f)
        return formation_energy

    def formation_energy_divacancies(self, nth = 1, index = 0, relax_params = {}):
        atoms = self.atoms.copy()
        atoms.calc = self.calc
        atom_energy = self.atom_energy
        Morph(atoms).create_divacancies(index1= index, nth = nth)
        relax(atoms, **relax_params)
        formation_energy = atoms.get_potential_energy() - atom_energy * len(atoms)

        dump_xyz('MaterialProperties.xyz', atoms)     
        with open('MaterialProperties.out', 'a') as f:
            f.write(f' {self.formula:<7}{nth}th Formation_Energy_Divacancies: {formation_energy:.4f} eV\n')
        return formation_energy

    def migration_energy_vacancy(self, index0 = 0, index1 = 1, fmax = 0.02, steps = 500):
        atoms = self.atoms.copy() 
        symbol = atoms[index0].symbol
        atoms[index1].symbol = symbol
        initial = atoms.copy()
        del initial[index0]
        final = atoms.copy()
        del final[index1]

        initial.calc = self.calc
        relax(initial)
        relax_cell = initial.get_cell()
        final.set_cell(relax_cell)
        final.calc = self.calc
        relax(final, method='ucf')

        images = [initial] + [initial.copy() for _ in range(11)] + [final]
        for i in images:
            i.calc = self.calc
        neb = NEB(images,climb=True, allow_shared_calculator=True)
        neb.interpolate()
        
        optimizer = FIRE(neb)
        optimizer.run(fmax=fmax, steps=steps)
        energies = [image.get_potential_energy() for image in images]
        energies = np.array(energies)
        migration_energy = max(energies) - min(energies)
        energies -= min(energies)
        for image in images:
            dump_xyz('MaterialProperties.xyz', image)  
        with open('MaterialProperties.out', 'a') as f:
            print(f' {self.formula:<7}Migration_Energy_({symbol}-Vacancy): {migration_energy:.4f} eV', file=f)
        plt.plot(np.linspace(0, 1, len(energies)), energies, marker='o', label=f'{self.formula}')  
        plt.legend()
        plt.savefig(f'{self.formula}_migration_{symbol}_vacancy.png')
        plt.close()
        return energies
    
    def formation_energy_sia(self, vector = (1, 0, 0), index = 0, relax_params = {}):
        atoms = self.atoms.copy()
        atoms.calc = self.calc
        atom_energy = self.atom_energy
        Morph(atoms).create_self_interstitial_atom(vector, index = index)
        relax(atoms, **relax_params)
        formation_energy = atoms.get_potential_energy() - atom_energy * len(atoms)

        dump_xyz('MaterialProperties.xyz', atoms)
        with open('MaterialProperties.out', 'a') as f:
            print(f' {self.formula:<7}{vector} Formation_Energy_Sia: {formation_energy:.4} eV', file=f)
        return formation_energy
    
    def migration_energy_sia(self, vector1, vector2, fmax=0.02, steps=500):
        atoms = self.atoms.copy() 
        initial = atoms.copy()
        final = atoms.copy()
        index = get_nth_nearest_neighbor_index(initial, 0, 1)
        Morph(initial).create_self_interstitial_atom(vector1, index = 0)
        Morph(final).create_self_interstitial_atom(vector2, index = index)

        initial.calc = self.calc
        relax(initial)
        relax_cell = initial.get_cell()
        final.set_cell(relax_cell)
        final.calc = self.calc
        relax(final, method='ucf')

        images = [initial] + [initial.copy() for i in range(11)] + [final]
        for i in images:
            i.calc = self.calc
        neb = NEB(images, allow_shared_calculator=True)
        neb.interpolate()
        
        optimizer = FIRE(neb)
        optimizer.run(fmax=fmax, steps=steps)
        energies = [image.get_potential_energy() for image in images]
        energies = np.array(energies)
        migration_energy = max(energies) - min(energies)
        energies -= min(energies)
        for image in images:
            dump_xyz('MaterialProperties.xyz', image)  
        with open('MaterialProperties.out', 'a') as f:
            print(f'{self.formula:^4}   Migration_Energy_SIA: {migration_energy:.4f} eV', file=f)
        plt.plot(np.linspace(0, 1, len(energies)), energies, marker='o')  
        plt.savefig(f'migration_SIA.png')
        plt.close()
        return energies
    
    def formation_energy_interstitial_atom(self, symbol, fractional_position, config_type, new_atom_e = 0, relax_params = {}):
        atoms = self.atoms.copy() 
        atoms.calc = self.calc
        atoms_energy = atoms.get_potential_energy()
        if new_atom_e == 0:
            new_atom_e = self.atom_energy
        position = np.dot(fractional_position, self.atoms.get_cell())
        insert_atom = Atom(symbol, position)
        atoms.append(insert_atom)
        atoms.calc = self.calc
        relax(atoms, **relax_params)
        formation_energy = atoms.get_potential_energy() - atoms_energy - new_atom_e

        dump_xyz('MaterialProperties.xyz', atoms)
        with open('MaterialProperties.out', 'a') as f:
            print(f' {self.formula:<7}{config_type} Formation_Energy: {formation_energy:.4f} eV', file=f)
        return formation_energy
    
    def migration_energy_interstitial(self, symbols, fractional_position, config_type, fmax = 0.02, steps = 500):
        atoms = self.atoms.copy()
        position0 = np.dot(fractional_position[0], self.atoms.get_cell())
        position1 = np.dot(fractional_position[1], self.atoms.get_cell())
        insert_atom1 = Atom(symbols[0], position0)
        insert_atom2 = Atom(symbols[1], position1)
        initial = atoms.copy()
        initial.append(insert_atom1)
        final = atoms.copy()
        final.append(insert_atom2)

        initial.calc = self.calc
        relax(initial)
        relax_cell = initial.get_cell()
        final.set_cell(relax_cell)
        final.calc = self.calc
        relax(final, method='ucf')

        images = [initial] + [initial.copy() for i in range(11)] + [final]
        for i in images:
            i.calc = self.calc
        neb = NEB(images, allow_shared_calculator=True)
        neb.interpolate()
        
        optimizer = FIRE(neb)
        optimizer.run(fmax=fmax, steps=steps)
        energies = [image.get_potential_energy() for image in images]
        energies = np.array(energies)
        migration_energy = max(energies) - min(energies)
        energies -= min(energies)
        for image in images:
            dump_xyz('MaterialProperties.xyz', image)  
        with open('MaterialProperties.out', 'a') as f:
            print(f' {self.formula:<7}   Migration_Energy_{config_type}: {migration_energy:.4f} eV', file=f)
        plt.plot(np.linspace(0, 1, len(energies)), energies, marker='o', label=f'{config_type}')  
        plt.legend()
        plt.savefig(f'migration_{config_type}.png')
        plt.close()
        return energies
    
    def stacking_fault(self, a, b, miller, distance):
        '''
        ---------------------------------------------------------------------------------------------------
        For FCC-Al                  |   For BCC-Nb
        surf. I        surf. II     |   surf. I              surf. II            surf. III
        (-1, 1,  0)    (-1, 1,  0)  |   (-1, 1, -2)          ( 1,-1,  0)         ( 1, -4/5, 1/5)
        ( 1, 1, -2)    ( 1, 1,  0)  |   (-1, 1,  1) <111>    ( 1, 1, -1) <111>   ( 1,  1,    -1)  <111>
        ( 1, 1,  1)    ( 0, 0,  1)  |   ( 1, 1,  0) {110}    ( 1, 1,  2) {112}   ( 1,  2,     3)  {123}
                                    
        ---------------------------------------------------------------------------------------------------
        For HCP-Ti 
        for basal {0001} Normal--[0001]          for prism {10-10} Normal--[10-10]              
        uvws =  [[-2,1,1,0],  [[-1, 0,0]          uvws =  [[-1,2,-1, 0],    [[0,1,0] 
                 [0,-1,1,0],   [-1,-2,0]                   [0, 0, 0, 1],     [0,0,1] 
                 [0, 0,0,1]]   [ 0, 0,1]]                  [1, 0,-1, 0]]     [2,1,0]] 
             
        ---------------------------------------------------------------------------------------------------                                                       
        for Pyramidal I narrow  {10-11} Normal: None    for Pyramidal I wide {10-11} Normal: None   
        uvws = [[-1, 2,-1, 0],  [0,  1,0]               uvws = [[-1,-1,  2, 3],     [[-1,-1,1]      
                [-1, 0, 1, 2],  [-2,-1,2]                       [-1, 0,  1, 2],      [-2,-1,2]      
                 {1,  0,-1, 1}]                                  {1,  0, -1, 1}]                 
                                                                                                                                                         
        for Pyramidal II {11-22} Normal       
        uvws =  [[-1,-1, 2,3],     [[-1,-1,1]
                [-1, 1, 0,0],      [-1, 1,0]
                {1,1,-2,  2}       [N,  N,N]]                                                                         
        ---------------------------------------------------------------------------------------------------
        Hexagonal Miller direction indices to Miller-Bravais indices and back:
        [-1,-1,2,3] = [-1,-1,1]
        [-1, 1,0,0] = [-1, 1,0]
        ---------------------------------------------------------------------------------------------------
        ''' 
        atoms = self.atoms.copy()
        atoms.calc = self.calc

        slab = cut(atoms, a, b, clength=None, origo=(0,0,0), nlayers = 18, extend=1.0, tolerance=0.01, maxatoms=None)
        rotate(slab, a,(1,0,0),b,(0,1,0), center=(0,0,0))
        slab.calc = self.calc
        relax(slab)
        slab.center(axis=2)
        slab.constraints = FixAtoms(indices=[atom.index for atom in slab if atom.position[2] < 1/2 * slab.cell[2][2]])
        box = slab.get_cell()
        S = (box[0][0] * box[1][1] - box[0][1] * box[1][0]) * 2

        shift_distance = np.linalg.norm(np.array(a)) * distance
        shift_indices = [atom.index for atom in slab if atom.position[2] > 1/2 * slab.cell[2][2]]
        slide_steps = shift_distance / 10
        
        energies = []
        for i in range(11):
            slab_shift = slab.copy()
            slab_shift.positions[shift_indices] += [slide_steps * i, 0,0]
            slab_shift.calc = self.calc
            relax(slab_shift, method='fixed_line')
            defects_energy = slab_shift.get_potential_energy() / S
            energies.append(defects_energy)
            dump_xyz('MaterialProperties.xyz', slab_shift)

        energies = [e - energies[0] for e in energies]
        with open('MaterialProperties.out', 'a') as f:
            print(f' {self.formula:<7}{miller} Stacking_Fault: {max(energies) * 1000:.4f} meV/A^2', file=f)

        plt.plot(np.linspace(0, 1, len(energies)), energies, marker='o', label=f'{self.formula}')  
        plt.legend()
        plt.savefig(f'{self.formula}_stacking_fault_{miller}.png')
        plt.close()
        return energies
    
    def bcc_metal_screw_dipole_move(self, fmax = 0.02, steps = 500):
        lc = self.lc
        symbols = []
        compositions = []
        for symbol, composition in re.findall('([A-Z][a-z]*)(\d*)', self.formula):
            symbols.append(symbol)
            compositions.append(int(composition) if composition else 1)
        if len(symbols) > 1:
            element_ratio = np.array(compositions) / sum(compositions)
            element_counts = np.ceil(element_ratio * 135).astype(int)
            symbols = np.repeat(symbols, element_counts)
            np.random.shuffle(symbols)
            sym = symbols[:135]
        else:
            sym = [symbols[0] for _ in range(135)]
        with tempfile.NamedTemporaryFile(mode='w') as f:
            f.write(unit_screw)
            initial_screw = read_xyz(f.name)
        for i in initial_screw:
            i.set_chemical_symbols(sym)
            i.calc = self.calc

        initial = initial_screw[0]
        final = initial_screw[1]

        unit_cell = initial.cell.copy()
        initial.set_cell(lc * unit_cell, scale_atoms=True)
        relax(initial, method='ucf')
        initial_cell = initial.cell.copy()
        final.set_cell(initial_cell, scale_atoms=True)
        relax(final, method='ucf')

        images = [initial] + [initial.copy() for i in range(15)] + [final]
        for i in images:
            i.calc = self.calc
        neb = NEB(images, allow_shared_calculator=True)
        neb.interpolate()    
        optimizer = FIRE(neb)
        optimizer.run(fmax=fmax, steps=steps)
        energies = [image.get_potential_energy() / 2  for image in images]
        energies = np.array(energies)
        energies -= min(energies)
        for image in images:
            dump_xyz('MaterialProperties.xyz', image)  

        plt.plot(np.linspace(0, 1, len(energies)), energies, marker='o', label=f'{self.formula}')
        plt.legend()
        plt.savefig(f'{self.formula}_screw_dipole_move.png')
        plt.close()
        return energies

    def bcc_metal_screw_one_move(self, fmax = 0.02, steps = 500):
        lc = self.lc
        symbols = []
        compositions = []
        for symbol, composition in re.findall('([A-Z][a-z]*)(\d*)', self.formula):
            symbols.append(symbol)
            compositions.append(int(composition) if composition else 1)
        if len(symbols) > 1:
            element_ratio = np.array(compositions) / sum(compositions)
            element_counts = np.ceil(element_ratio * 135).astype(int)
            symbols = np.repeat(symbols, element_counts)
            np.random.shuffle(symbols)
            sym = symbols[:135]
        else:
            sym = [symbols[0] for _ in range(135)]
        with tempfile.NamedTemporaryFile(mode='w') as f:
            f.write(unit_screw)
            initial_screw = read_xyz(f.name)
        for i in initial_screw:
            i.set_chemical_symbols(sym)
            i.calc = self.calc

        initial = initial_screw[0]
        final = initial_screw[2]

        unit_cell = initial.cell.copy()
        initial.set_cell(lc * unit_cell, scale_atoms=True)
        relax(initial, method='ucf')
        initial_cell = initial.cell.copy()
        final.set_cell(initial_cell, scale_atoms=True)
        relax(final, method='ucf')

        images = [initial] + [initial.copy() for i in range(15)] + [final]
        for i in images:
            i.calc = self.calc
        neb = NEB(images, allow_shared_calculator=True)
        neb.interpolate()    
        optimizer = FIRE(neb)
        optimizer.run(fmax=fmax, steps=steps)
        energies = [image.get_potential_energy() for image in images]
        energies = np.array(energies)
        energies -= min(energies)
        for image in images:
            dump_xyz('MaterialProperties.xyz', image)  

        plt.plot(np.linspace(0, 1, len(energies)), energies, marker='o', label=f'{self.formula}')
        plt.legend()
        plt.savefig(f'{self.formula}_screw_one_move.png')
        plt.close()
        return energies